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ABSTRACT 

In a systematic study, we compare the density statistics in high-resolution numerical experiments of 
supersonic isothermal turbulence, driven by the usually adopted solenoidal (divergence-free) forcing 
and by compressive (curl-free) forcing. We find that for the same rms Mach number, compressive 
forcing produces much stronger density enhancements and larger voids compared to solenoidal forcing. 
Consequently, the Fourier spectra of density fluctuations are significantly steeper. This result is 
confirmed using the A-variance analysis, which yields power-law exponents (3 ~ 3.4 for compressive 
forcing and [3 ~ 2.8 for solenoidal forcing. We obtain fractal dimension estimates from the density 
spectra and A-variance scaling, and by using the box counting, mass size and perimeter area methods 
applied to the volumetric data, projections and slices of our turbulent density fields. Our results 
suggest that compressive forcing yields fractal dimensions significantly smaller compared to solenoidal 
forcing. However, the actual values depend sensitively on the adopted method, with the most reliable 
estimates based on the A-variance, or equivalently, on Fourier spectra. Using these methods, we 
obtain D ~ 2.3 for compressive and D ^ 2.6 for solenoidal forcing, which is within the range of 
fractal dimension estimates inferred from observations {D ^ 2.0 . . . 2.7). The velocity dispersion to 
size relations for both solenoidal and compressive forcings obtained from velocity spectra follow a 
power law with exponents in the range 0.4 . . . 0.5, in good agreement with previous studies. 
Subject headings: hydrodynamics — ISM: clouds — ISM: kinematics and dynamics — ISM: structure 
— methods: numerical — turbulence 
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1. INTRODUCTION 

Observations provide velocity dispersion to size rela- 
tions for various molecular clouds (MCs), which doc- 
ument the existence of supersonic ra ndom motions 
on scales larger than ^ 0.1 pc (e.g., Larson 
Myers 1983; Pcrault ct al. 1986; Solomon et al. 1987; 
Falgarone et ail Il99^ : iHever fc BruntI 12004ft . These 



moti ons are associated with compressible turbulence 
(e.g.. lElmegreen fc Scaloll2004H Scal o fc ElmegreCTl l2004 
IMac Low fc KlessenI l20oir in the interstellar medium 
(iFerrierd |2001D and exhibit a single turbulent cas- 
cade or spat i ally s eparated coexisting inertial ranges 
(|Passot et al.l Il988f ) sin iilar to the k inetic energy cas- 
cade of incompressible iKolmogorovl (|1941h turbulence. 
However, there are various physical processes (e.g., self- 
gravity, magnetic fields, nonequilibrium chemistry) and 
especially the cor npressibility of the gas, that alter the 
scafing laws (e.g.. lFlec3ll996D an d statistics (e.g., inter- 
mittency corrections measured bv lHilv-Blant et a"l]|2008f) 
established for incompressible turbulence. 

The physical origin and characteristics of the turbu- 
lent fiuctuations are still a matter of debate. To advance 
on the question of how turbulence isdriven in the inter- 
stellar medium, we present results of high-resolution nu- 
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merical experiments of supersonic isothermal turbulence 
comparing two distinct and extreme ways of driving the 
turbulence in a systematic study: 1) solenoidal forcing 
(divergence-free or rotational forcing), and 2) compres- 
sive forcing (curl- free or dilatational forcing). 

Various numerical and analytical studies have 
provided important insight into the statistics of 
super s onic isothermal t urbulence (e.g.. iPorter et al.l 



1992t Ivizauez-Semadenil Il99l iPadoan et all 
Passot fc Vazquez-Semadenil 119981: 'Stone et al.l 



Mac Low! 1999: Klessen 2000; Ostriker et all 



Boldvrev et al.ii2002i: iLi et al. 2003; Pado an et al 112004 



Jappsen et al.l 2005 : B allcs tcros-Parc des et al. I 1200 



Kritsuk et aTri2007l : (Lemaster 'fc Stonell2008f) . Most of 



these studies use purely solenoidal or weakly compres- 
sive kinetic energy injection mechanisms (forcing) to ex- 
cite turbulent motions. In the present study, we aim 
at comparing the usual case of solenoidal (divergence- 
free) forcing with the case of fully compressive (curl- 
free) forcing. The actual way of turbulence production 
in real MCs is expected to be far more complex com- 
pared to what we can model with the present simula- 
tions, probably consisting of a convolution of various 
agents producing turbulence, and mixtures of s olenoidal 
and compressive modes ( e.g., lElmegreen fc Scalo. ,200^ 
IMac Low fc KlessenI l200'3) . Here, we systematically in- 
vestigate the extreme cases of purely solenoidal versus 
purely compressive energy injection. 

Analyzing the density correlation statistics and frac- 
tal structure obtained in our hydrodynamic simulations, 
we show that compressive forcing leads to significantly 
steeper density fluctuation spectra and consequently to 
fractal dimensions of the turbulent gas structures, that 
are significantly smaller compared to the usually adopted 
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solenoidal forcing. We use Fourier analysis, A-variance 
analysis, structure functions, the fractal mass size, box 
counting, and perimeter area methods to obtain fractal 
dimension estimates. We apply the A-variance analy- 
sis to both our 3-dimensional data and to 2-dimensional 
projections, and the perimeter area method to projec- 
tions and slices through the turbulent density structures 
supporting the result of a significantly smaller fractal di- 
mension for compressive forcing compared to solenoidal 
forcing. Although compressive forcing yields signifi- 
cantly smaller fractal dimensions than solenoidal forc- 
ing, the estimated fractal dimensions are in the range 
2.0... 2. 7 consistent with observational estimates (e.g., 
lElmegreen fc Falgaronj|1996HSanchez et al.ll2007[ ) 

We explain our numerical method, construction of 
solenoidal and compressive forcing fields and fractal anal- 
ysis techniques in Section [51 In Section [31 we show 
that our results are consistent with previous studies us- 
ing solenoidal forcing, whereas compressive forcing yields 
much stronger density contrasts and consequently leads 
to significantly smaller fractal dimensions. In Section [H 
we summarize our conclusions. 



2. SIMULATIONS AND METHODS 

Th e piecewise parabolic method (Colell a fc WoodwardI 
|1984|) implementati o n of the astrophysi cal code FLASH3 
([Frvxell et al] 120001 : iDubev et all 120081) was used to in- 
tegrate the hydrodynamic equations on periodic uniform 
grids with 256^, 512^ and 1024^ grid points. Density p, 
velocity v and pressure P are related through the equa- 
tions 

' - (pv)=0 (1) 



„^ , V-(pv)=0 
ot 

^ + (v.V)v = -ivP + f . (2) 
ot p 

Note that an energy equation is not needed, because we 
model isothermal gas. The pressure is simply given by 
P = cj.p with the constant sound speed Cg. Isother- 
mality is a very crude, but reasonable first approxima- 
tio n for modeling the thermodynarn ic behavior of MCs 
IWolfire et al.l 119951: iPavlovski et al.r 2006l. Due to the 
isothermal approximation, the hydrodynamic equations 
are scale-free, and we can solve them for a chosen den- 
sity scale po = 1 (mean density), sound speed Cg = 1 and 
domain size L — 1. The only remaining free parame- 
ter therefore is the dimensionless rms Mach number A4, 
which can be varied. It is important to note that the 
forcing term f used to drive turbulent motions appear- 
ing as source term in equation ([2]) can also be varied. In 
the present study, we vary the forcing term, investigat- 
ing the difference between purely solenoidal and purely 
compressive kinetic energy injection, while keeping the 
rms Mach number fixed. 

Equations ^ and ([2]) have been solved 
merically with periodic boundary conditions 
ing an isothermal equation of state in the 
text of MC dyn amics in various studies 
Padoan et al. 1997 1: iPassot fc Vazauez-Semadenil 



Stone et al.lil9 98: Mac Low et al."1998l; iMac Lowl 
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2002; Li ct al. 2003; Padoa n et al.ll2004bl : IJappsen et al 



2005; .Ballesteros-Paredes et al.l l2006t 
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comparing turbulence statistics obtained in our study 
with the results of these studies. In particular, we want 
to check the influence of different forcings. Therefore, we 
concentrate on two extreme cases: 1) the usually adopted 
solenoidal forcing (divergence-free forcing) and 2) fully 
compressive forcing (curl-free forcing). 

2.1. Forcing Module 

Turbulent fluctuations have to be excited and main- 
tained in order to study stationary turbulence statistics 
in detail. If not constantly driven by a random force 
field, turbulent motions damp due to dissipation. In most 
studies, the force field is constructed in Fourier spa ce by a 
3-di mensional stochastic procedure (e.g., lOu binski et aTl 
1995; Mac Low et"aIlll998l:IStone et al.lll998f ). which gen- 
erates a random vector field f after Fourier transfor- 
mation back into physical space. This field will on 
average contain 2/3 of its energy in solenoidal modes 
(transversal modes) and 1/3 in compressive modes (lon- 
gitudinal modes), because in 3-dimensional space, waves 
have two spatial directions for the transversal part, 
whereas the longitudina l part has only one (see, e.g., 
lElmegreen fc Scaloll2004[ ). In order to obtain a purely 
solenoidal, or a purely compressive forcing field f, a 
Helmholtz decomposition can be made by applying the 
projection operator vjj in Fourier space (/c-space) 



cr^ + (1 - C)vl = cs^J + (1 - 20 



W 



(3) 



We aim at 



prior to the inverse Fourier transformation into real 
space. By setting the parameter C G [0, 1] one can adjust 
the mixture of solenoidal and compressive modes. If we 
set C = 1, V'lj projects only the solenoidal component, 
whereas only the compressive component is obtained by 
setting C = 0- 

The forcing term f is typically either modeled as 
a spatially static pattern with ti me-dependent ampli- 
tude (following th e recipes, e.g., bv lMac Low et al]|1998l : 
IStone et al.l I1998D or by usi ng an Ornstcin-Uhlenbeck 
(OU) process fe.g.. lEswaran fc Pope 1988; Schmidt ct al] 
l2006f ). which modulates the pattern smoothly in space 
and time on a well-defined autocorrelation timescale T 
resulting in a constant energy input rate. We follow 
the usual approach and set the autocorrelation timescale 
equal to the dynamical timescale T — L/{2V), where 
L is the size of the computational domain, V = CsA4 
and M ~ 5.5 is the rms Mach number in all runs. 
Therefore, T is the time for the most energetic fluctu- 
ations (at fc = 2 in Fourier space, which corresponds to 
L/2) to cross half of the box. It is furthermore equal to 
the decay time consta nt of the turbulence (Stone ct af] 
119981: IMac Low! [l999l ). The forcing amplitude follows 
a parabolic power spectrum only containing power on 
the largest scales in a small interval of wavenumbers 
1 < fc < 3 peaking at A; = 2. The influence of varying the 
sc ale of ener g y inp u t has been inv estigated fo r instance 
bvlMac Low! (Il999f). [Klessen et al.l (120001 Heitsch et al.l 
(|2001h and IVazquez-Semadeni et all (|2003i ). Here, we 
only consider the usually applied large-scale stochastic 
forcing. This way of forcing models the kinetic energy 
input from larger scale turbulent fluctuations breaking 
up into smaller structures and feeding kinetic energy to 
smaller scales. 
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We checked that our results are not sensitive to the 
particular method for generating turbulent motions, i.e., 
by using an almost static pattern (using a very large au- 
tocorrelation time in the OU process), and by using a 
band spectrum instead of a parabolic Fourier spectrum 
for forcing. Variations in the spectral form of the large 
scale forcing did not significantly change the results ob- 
tained in the present study. However, changing the mix- 
ture of modes from a purely solenoidal to a purely com- 
pressive forcing always yielded significant differences. 

2.2. Initial Conditions and Post Processing 

Starting from a uniform density distribution and zero 
velocity, the forcing excites turbulent motions. The forc- 
ing amplitude is adjusted to excite turbulence with rms 
Mach number Ai ~ 5.5. We use Ai as the control pa- 
rameter, because this dimensionless number is often ex- 
pected to solely determine physical properties of scale- 
invariant turbulent flows. The purpose of the present 
study was to determine the effect of varying the forcing 
from purely solenoidal to purely compressive, so we keep 
the rms Mach number fixed besides all other parame- 
ters and varied only the forcing between solenoidal and 
compressive. 

Equations ([T]) and ([2]) were evolved for ten dynamical 
timescales T, which allows us to study a large sample of 
statistically stationary realizations of the turbulent flow. 
We wait for two dynamical timescales before averaging all 
statistical measures in the time interval 2 < t/T < 10. 
Since we have produced snapshots every 0.1 T, the re- 
sulting statistical sample consists of 81 realizations of 
the turbulent field. The averaging procedure is impor- 
tant to derive meaningful statistics, becaus e all quantities 
are su bject to statistical fluctuations (e.g.. lKritsuk et all 
I2OOI . The averaging procedure furthermore provides a 
handle on the Icr temporal fluctuations between differ- 
ent realizations. Unless otherwise stated, the Icr tempo- 
ral fluctuations are indicated as error bars in the results 
section. 

2.3. Box Counting Method 

We analyzed fractal structures in our simulation data 
using the box counting method. In the first step, the frac- 
tal structure is defined by marking all cells belonging to 
the fractal set, if they are above a certain density thresh- 
old, whereas all cells below that threshold are marked as 
not belonging to the fractal structure. In the second step, 
the structure as defined above was scanned by applying 
a box (mask) of size I and counting how often the struc- 
ture is covered by that box. This procedure was repeated 
varying the size of the box resulting in a set of counts iV^ 
and box sizes k. A plot of log(A'i) against log(/i) often 
reveals a scaling range over which the points fall close 
to a straight line with the box cou nting dimension Dh as 
the negative slope of that line (e.g.. [Mandelbrot &: Frame! 
l2002tlPeitgen et al.|[200i ). This implies a power-law scal- 
ing N{1) oc in the scaling range. 

Setting the density threshold pth for defining the frac- 
tal structure is obviously a critical choice. Using pth = 
naturally results in Db — 3, whereas setting pth = Pmax 
leads to Dt, — 0. We computed the box counting di- 
mension for different density thresholds and discuss its 
influence on the results. 



2.4. Mass Size Method 

The fractal mass size dimension was obtained by com- 
puting the mass contained inside concentric boxes with 
increasing box size / centered on the densest cells of the 
data set and averaging over cells with p > Pmax/2, fo l- 
lowing the method described by iKritsuk et al.l (|2007f) . 
This yields a set of masses Mi and box sizes k. A plot 
of log(Mi) against log(Zi) often reveals a scaling range 
over which the points fall close to a straight line with the 
mass size dimension Dm as the slope of that lin e (e.g., 
[Mandelbrot fc Framd l2002t iPeitgen et aP l2004( l. This 
implies a power-law scaling M{1) (x Z^™ in the scaling 
range. 

Power-law relations of the form M{1) cx Z° should be 
considered with caution in the context of fractals, be- 
cause such relations often occur in physics and do not 
necessarily imply that a is a fractal dimension, i.e., for a 
3-dimensional d ensity d istribution, a > 3 can o ccur (see 
lMandelbrotl[l983: Elmegreen fc Falgaronelll996l ). 

2.5. Perimeter Area Method 

We also applied the perimeter area method, which is 
frequently applied for measuring the fractal dimension 
of interstellar gas clouds. The boundary curves and ar- 
eas of coherent structures with equal density were identi- 
fied in both, 2-dimensional projections and 2-dimensional 
slices through the computational domain. Varying the 
density threshold yields a set of structures with perime- 
ters Vi and areas Ai- Fitting a power law of the form 
V oc A^''^^ (log- log plot as for the box counting and 
mass size dimensions) yields the perimeter area dimen- 
sion Dp. Structures with very smooth boundary curves 
exhibit Dp = 1, whereas for structures with totally con- 
voluted perimeters, V grows linearly with the area occu- 
pied by the structure resulting in Dp = 2. 

3. RESULTS AND DISCUSSION 

3.1. Time Evolution 

Figure [T] compares projections (top panels) and slices 
(bottom panels) of the density field in the x-y-plane from 
a randomly picked snapshot (t = 5 T) for solenoidal 
versus compressive forcing as an example of the typi- 
cal density structure in the state of statistically station- 
ary supersonic turbulent fiow. This regime was safely 
reached after 2 dynamical times T, which is demon- 
strated in Figure [2l The rms Mach number has settled 
to ~ 5.5 for both solenoidal and compressive forcing, 
and for numerical resolutions of 256'^, 512'^ and 1024'^ 
grid points after 2 dynamical times. Not only the ve- 
locity statistics has converged to a stationary state, but 
also the density statistics, which is shown in terms of 
minimum and maximum densities in the top panel of 
Figure [21 Obviously, compressive forcing produces larger 
density contrasts, which results in higher density peaks 
and larger voids. Although both cases exhibit the same 
rms Mach number, the solenoidally driven case gives a 
much smoother density distribution with smaller disper- 
sion. In both cases, the maximum de nsity is subject 
to strong intermittent fiu ctuations fe.g. jFalgarone et al.l 
1 1994 IKritsuk et al.ll200"7b leading to temporal variations 
in the maximum density of order one magnitude. 

3.2. Fourier Spectrum Functions 
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and density fluctuation Fourier spectra 

P{k)dk= / (/r^)(p^^o)*47rfc2dfc. 



We begin the analysis of spatial correlations by show- 
ing 1-dimensional Fourier spectra derived for solenoidal 
and compressive forcing. Figure [3] presents a comparison 
of velocity Fourier spectra 

E{k) dk^- / V • 9* iirk^dk (4) 



(5) 



These were computed with the standard method (e.g., 
iFrischI |1995[ ). i.e., by integration over spherical shells 
in Fourier space using the Fourier transformed velocity 
components Vi{k) and the Fourier transformed density 
fluctuations p — po(k), where po denotes the mean den- 
sity. Velocity Fourier spectra E(k) are ty pically used 
to distinguish between iKolmogorovl ()194ll ) turbulence, 
E{k) oc k~^/^ and Burgers turbulence, E{k) oc in 
the inertial range. For highly compressible isothermal 
supersonic turbulent flow, it has been shown that the in- 
ertial ran ge scaling is c l oser to Burgers turbulence. For 



instance, 



iKritsuk et all (l200l find Eik) cx fc-^'^^ and 

iSchmidt et all (|200 9^ measured E{k) oc k In the 

present study, we obtain E{k) oc k~^'^'^ for compressive 
forcing, slightly steeper than E{k) oc fc~^ ®^ for solenoidal 
forcing. Taking the error bars (temporal variations) into 
account, this represents just marginal difference between 
both forcing schemes, and our estim ates with i n the iner- 
tial range are consistent with Kritsu k~et al.l ()2007f ) and 
ISchmidt et all (p00 9). Table 1 summarizes all results 
obtained for solenoidal and compressive forcing provid- 
ing the formal least-squares fitting errors, which were 
obtained by taking into account the la temporal fluc- 
tuations. Table 1 furthermore provides a summary of 
defining relations for the following fractal dimension es- 
timates. 

Note that simila r to IKritsuk et al.l (|2007t l and 
ISchmidt et all (l2009l) . we define our inertial range in 
a very small range of wavenumbers 5 < /c < 15 be- 
cause even at numerical resolutions of 1024'^ grid points. 



the i nertial range is very small (see, e.g. . .Klein et al 
20071 ) due to the bottleneck effect fe.g.. iDobler et al 



gOOa lHaugen fc Br andenburgll2004l : [Schmidt et al.ll2006l : 

Kritsuk et al. 2007), which affects the Fourier spectrum 
in the dissipation range. For our simulations, we esti- 
mate that wavenumbers fc > 20 may be affected by the 
bottleneck. The infiuence of the bottleneck effect is less 
pronounced in physical space, which allows for a some- 
what larger fitting range for scaling estimates obtained 
in physical space (simila r to, e.g.. iKowal fc LazarianI 
[20071: kritsuk et all 120071 : ISchmidt et all [20f)9h for the 
A- variance analysis in Section 15?^ 

The density fiuctuation spectra in Figure [3| show con- 
siderable differences in their amplitude and inertial range 
scaling for solenoidal and compressive forcing. First, 
we discu ss the difference in th e amplitudes. As dis- 
cussed in lFederrath et all ()2008f ). the standard deviation 
of the density probability distribution function (PDF) is 
~ three times larger for compressive forcing compared to 
solenoidal forcing. This result is recovered in the present 
analysis by integrating the density fluctuation spectra 



where n — 1024'^ is the number of grid points, which 
yields the standard deviation a p. For solenoidal forcing, 
we compute J P{k)dk ~ 1.89, whereas for compressive 
forcing, J P{k)dk ~ 5.93, in very good agreement with 
the standard devi ations go^l.89 and ao ~5.86 obtained 
from the PDFs in lFederrath et al.l (|2008h . 

Second, the inertial range scaling differs significantly 
for the two cases. For solenoidal forcing, P{k) oc fc"" ''^ 
and for compressive forcing, P{k) oc k~^ '^^. Our iner- 
tial range scaling inferred for solenoidal forcing is con- 
sistent with the weakly ma gnetized sup e r-Alfv enic su- 
personic MHD models f romlKowal et al.l (l2007l Tab. 2) 
using solenoidal forcing. iKowal et al.l ( 2007f l find P(fc) oc 
^-o.6±o.2 jtqj, tjjeij- model B.lP.Ol with very weak mag- 
netic field, slightly larger rms Mach number [AA ^--^7.1), 
and resolution'* of 256'^ explaining the small differences 
comparing their result with ours. It is furthermore in 
agree ment wi th the purely hydrodynamic estimates by 
Krits uk et al.l (2006) with resolutions up to 2048'^ using 
adaptive mesh refinement, who obtain P{k) oc fc-0-8 - 0.9^ 

Note that in general, the power-law exponents a of 1- 
dimensional Fourier spectra are related to the power-law 
exponents (3 of the corresponding 3-dimensional Fourier 
power spectra thro ugh f3 = a + 2. As discussed by 
IStutzki et aT] (|1998f). the power-law scaling of the density 
spectrum is furthermore related to the fractal drift expo- 
nent H (Hurst exponent). Considering a 1-dimensional 
density power spectrum P{k) oc fc~" leads to a Hurst 
exponent oi H = a/2 = (/? - 2)/2. The Hurst expo- 
nent is related to the fractal box coverage dimension 
Z3b = E + 1 — H with the dimensionality = 2 for 
the box coverage of a fractal surface embedded in 3- 
dimensional space ijStutzki et al.|[l998[ ). Consequently, 
we obtain H ^ 0.39 and Di, ~ 2.61 for solenoidal forc- 
ing, and H^Q.72 and Z?fc~2.28 for compressive forcing. 
Thus, the fractal Hurst exponent is significantly smaller 
for compressive forcing. The estimates for the Hurst ex- 
ponents and the corresponding relations with box count- 
ing and perimeter area dimensions are summarized in 
Table 1. 

3.3. A- Variance Analysis 

In this section, we present results of the A-variance 
method, which provides an independent measure of the 
scaling of the density Fourier spectra. IStutzki et al.l 
(1998) introduced a wavelet transformation called A- 
variance for measuring the exponent of the Fourier spec- 
trum. As complementary analysis, we performed the 
A-variance with the t ool developed and provided by 
[Ossenkopf etin (|2008D^t_ is an im proved vers ion of 
the original A-variance (jStutzki et al.f l998: Bcnsch et al] 
l200lf l. The A-variance measures the amount of structure 
on a given length scale I, filtering the data set <i>(x) by ap- 
plying an up- down- function Qi (typically a French-hat 
or Mexican-hat filter) of size I and computing the vari- 
ance of the filtered data set. The A-variance is defined 
as 



^a(0 



$(x) * Q (x; 



(7) 



P{k)dk, 



(6) 



where the average is computed over all valid data points 
at positions x, and the operator * stands for the con- 

* Th e res olution dependence for our simulations is discussed in 
Section 13.31 
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volution. The data set can have arbitrary dimension- 
ality. In the present study, we apply the A-variance 
to both, 2-dimensional (projections) and 3-dimensional 
data sets. We checked the influence of varying the fil- 
ter function from French-hat to Mexican-hat, as well as 
varying the diameter ratio of the filter, which yielded 
no significant differences. The choice of the filter func- 
tion and t he best choice f o r its diameter ratio is dis- 
cussed by lOssenkopf et al.l (|2008D . Here, we use the 
original French-hat filter with a d iameter ratio of 3.0 
as it has been used b efore (e.g., ISt utzki et al.l 1X998 



Mac Low fc Ossenkopfl |200(J : [Osscnko pf et al 
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Ossenkopf fc Mac Lo^ 120021: lOssenkopf et all l2006f ). 
Note that originallv. IStutzki et al.l (|1998D applied the 
A-variance to 2-dimensional integrated maps for com- 
parison with observations. Although we have access to 
the 3-dimensional data from our simulations, we never- 
theless computed column density maps and applied the 
A-variance to both the 2-dimensional and 3-dimensional 
data to determine the effect of projection for applying 
the A-variance. Prior to the 3-dimensional analysis, we 
resampled the density data cubes with 1024'^ grid points 
to a resolution of 256^ due to performance issues of the 
A-variance, which is not (yet) a parallelized tool. The 
resampling to lower resolution is not expected to cause 
deviations in the inertial range scaling. As long as the 
simulation itself had enough spatial resolution to resolve 
the inertial range scaling, the resampling to lower reso- 
lution only affects the dissipation range leaving density 
spectra alm ost up to the new Nyquist frequency (e.g., 
iPadoan eF al. 2006) . We explicitly show in the bottom 
panel of Figure [4] that the resampling indeed did not 
affect the inertial range. Only the compressive forcing 
case is shown but the resampling for the solenoidal case 
exhibits similar behavior. 

The upper panel of Figure |4] shows the influence of 
varying the numerical resolution of the simulation. The 
inertial range scaling depends on the numerical resolu- 
tion. A resolution of 256^ grid points seems insufficient 
to resolve the exact inertial range scaling, although the 
15% difference compared to the 1024'^ simulation is of the 
order of the temporal fluctuations, whereas the difference 
between solenoidal and compressive forcing (Fig. [3|) is 
significant. At resolutions of 512^ and 1024'^ grid points, 
the best-fit power-law scaling agrees quite well, indicat- 
ing almost numerical convergence. A similar conclusion 
ca n be drawn from the density Fourier spectra presented 
by iKritsuk et al.l (|2006l Fig. 4) computed for solenoidal 
forcing. 

The results of the A-variance are presented in Fig- 
ured In the top panel, we show the A-variance applied 
to the 3-dimensional data resampled to 256^ grid cells, 
whereas the bottom panel presents the A-variance ap- 
plied to projections averaged along all three spatial axes. 
The variation for different projections is almost negligible 
compared to the temporal fluctuations. We nevertheless 
added the variation due to projection along the three dif- 
ferent axes to the \a error bars due to temporal fluctu- 
ations. Following IStutzki et al.l (|1998[ ). the values of the 
best-fit power-law exponents (3 are shown corresponding 
to the 3-dimensional Fourier spectra. Note that the slope 
a fitted to the A-variance is related to the slope of the 3- 
dimensional Fourier spectrum by /? = a+E, where E = 2 
for the projected data and E ~ 3 for the 3-dimensional 



data resu lting in the s a me p ower-law exponent /3. As 
shown bv IStutzki et al.l (|1998f ). the power law scaling of 
the Fourier spectrum should be preserved at projection 
to lower dimensions. I n agre ement with the results by 
iMac Low fc Ossenkopfl (2000), we find that the slopes 
/3-2.89 (3D) and /3~2.81 (2D projection) for solenoidal 
forcing, and f3 ~ 3.44 (3D) and /3 ~ 3.37 (2D projection) 
for compressive forcing are almost preserved during pro- 
jection (see Table 1). 

We can furthermore check whether the A-variance re- 
sults agree with the Fourier power spectra, since the 
A-variance is supposed to measure the power-law expo- 
nent of the Fourier spectrum. As shown in Figure [31 the 
density spectra exhibit power laws in the inertial range 
corresponding to 3-dimensional Fourier power-law expo- 
nents f3 = a + 2 = 2.78 for solenoidal and f3 = 3.44 for 
compressive forcing in consistency with the A-variance. 
Therefore, the A-variance confirms the results obtained 
by the density Fourier spectra, showing that compressive 
forcing yields a systematically steeper density spectrum 
compared to solenoidal forcing. 

3.4. Structure Functions 

Besides the Fourier spectra and the A-variance 
analyzed in the previous sections, structure functions 
are often used to measure spatial corr elations in tur- 
bulen t v elocity and de nsi ty fie lds (e.g., iBoldvrev et al.l 



lEsauivel fc 



2002[ IPadoan et al l 120031 : 
2005t IKritsuk et all 120071: iSchmidt et a 
Hilv-Blant et al.l 120081) . Here, we restrict our anal 



azariani 
120091 : 



ysis to the 2nd order structure functions of the density 
field for comparison with the Fourier spectra and 
A-variances. The 2nd order structure function of the 
density is defined as 



(8) 



SF2(0- |p(x)-p(x + l)r 



Figure [6] presents the 2nd order density structure func- 
tions for compressive and solenoidal forcing respectively. 
In the following, we draw the connection of these struc- 
ture functions to the power spectra and A-variance. One 
feature of the structure func tion is its relation t o the au- 
tocorrelation function A(l) (jStutzki et al.|[T998[) : 



SF2(0 = 2 [A(0) ~ Ail)] = 2 _ Ail)] 



(9) 



Since the autocorrelation function vanishes on large 
scales close to our periodic box size {1 ^ L), the 2nd or- 
der structure function of a variable can be used to mea- 
sure the standard deviation a of this variable because 
SF2(Z L) = 2(7^. In our case, we obtained the stan- 
dard deviations of the density Cp ^1.88 for solenoidal and 
(Tp ~ 5.95 for compressive forcing simply by evaluating 
cTp = [0.5SF2(Z = 0.5i)]i/2 from Figure El Note that 
for periodic boxes, the maximum length scale for mea- 
suring spatial correlations is half of the box size L. As for 
the power spectra, this is in good agree ment with the in- 
depen dent analysis of the density PDFs (jFederrath et al.l 
I20n8h . 

The best-fit values of power-law exponents SF2(Z) oc l^ 
of the structure functions in the inertial range are in- 
dicated in Figure |6| as well. Since SF2(?) oc l^ is the 
Fourier transform of the 1-dimcnsional Fourier spectrum 
P(fc) oc fc^", it follows that a = r/ + 1. For compressive 
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forcing, the power-law scahng is in agreement with the A- 
variance and Fourier spectra estimates. The correspond- 
ing value for the 3-dimensional density power spectrum 
derived from the structure function is/3 = a-|-2 = ?7-|-3~ 
3.47 (see Table 1). For solenoidal forcing on the other 
hand, the best-fit value derived from the structure func- 
tions is ~3.24 is in significant disagreement with the esti- 
mates from the Four ier spectra and A- variance (/3~2.8). 
iStutzki et al.l (|1998f ) provide the mathematical explana- 
tion for this. For a power-law scaling of Fourier spectrum 
functions with power-law exponent (3 < E (here i? = 3), 
i.e. a < 1, the £'-dimensional two-point correlation func- 
tion (structure functi on) does not necessa rily increase in 
a power-law fashion (jStutzki et al.l[l9 98L Appendix B). 
This limits the applicability of structure functions for 
estimating the inertial range scaling to density Fourier 
power-law exponents E < f3 < E + 2. 

3.5. Mass Size Method 

Figure [7] shows the results obtained by applying the 
mass size method as described in Section 12.41 In rough 
agreement with the results of the methods discussed so 
far, compressive forcing yields a smaller mass size ex- 
ponent Drn ^ 2.03 compared to solenoidal forcing with 
Dm ~ 2.11 in the inertial range (Table 1). Unlike the 
other methods, however, this difference is not significant. 
The large la error is a direct consequence of the strong 
temporal fluctuations of the maximum density seen in 
Figure [H Since the mass size relation M{1) is com- 
puted by growing concentric boxes centered on density 
peaks with p > Pinax/2, the mass is expected to fluctuate 
strongly, following the temporal fluctuations of the den- 
sity peaks. Our results are th erefore roughly consistent 
wit h the mass size analysis by iKowal fc Lazarianl ()2007f ) 
and lKritsuk et al.l (|2007[ ). 

3.6. Box Counting Method 

The results of the analysis using the box counting 
method as explained in Section 12.31 are presented in Fig- 
ure m In this case, the fractal structure was defined by 
the mean density po = 1 as threshold density for both 
solenoidal and compressive forcing. The analyzed struc- 
ture as defined by po is shown in Figure [9] We obtain 
the box counting dimensions Di, ~ 2.68 for solenoidal, 
and Df, ~ 2.51 for compressive forcing in the inertial 
range (Table 1). As discussed in Section [2.31 the box 
counting dimension depends on the threshold for defin- 
ing the fractal structure. We have checked its depen- 
dence on the threshold pth by varying r = log]^Q(/9th//3o)- 
The results obtained by computing the box counting di- 
mension for T = —1, 0, 1, 2 arc shown in Figure [TOl for 
solenoidal forcing (left panel) and compressive forcing 
(right panel). Note that r = corresponds to the mean 
density as defining threshold. As expected, the box di- 
mension strongly depends on pth- Significant differences 
between solenoidal and compressive forcing are obtained 
for different threshold densities. For thresholds t > 0, 
the box dimension is smaller for solenoidal forcing than 
for compressive forcing in contrast to r < 0. The latter 
is as a consequence of the much more space filling den- 
sity structure for the solenoidal case (see Fig. [1]). On 
the other hand, for r = 2 the solenoidal case yields 
filamentary structures with small fractal box dimension 
{Df,^0.72), while structures in the compressive case are 



coherent elongated and almost sheetlike structures with 
larger fractal dimension {Di, ~ 1.63). Obtaining abso- 
lute estimates for the fractal dimension using the box 
counting method in the present context is consequently 
impossible. However, differences among different data 
sets, e.g., solenoidal versus compressive forcing can be 
measured with the box counting method, if the same 
defining density threshold is used. 

3.7. Perimeter Area Method 

In this section, we show results of the perimeter 
area method described in Section 12.51 This method 
is often applied to observational data sets to infer 
the fractal dimension of interstellar clouds. Although 
we are aware of the fact that our simulations can 
only provide a very limited approximation to real 
interstellar gas, we nevertheless are convinced that 
comparison with observational data will provide phys- 
ical insight. The perimeters of interstellar gas clouds 
exhib it fractal di r nensions in the range Dp ^ 1.2 . . . 1.6 
Ce.g. lBeech'1987^ 'Bazell fc Deserti'igSB": i Dickman et al.l 



1990i Falgaronc et al. 1991 
Hetem & Lcoinc 1995; 



Voeelaa r et al 



1991 

Wakker 1994 

1999; iKim et al.l |2003; ,Lee ,2004 



fVogeL 



aar 



Westpfahl et al.. . . , . . , . . _ _ . 

Sanchez et al. r i2007D with most of the studies finding 



D„ 



1.3... 1.4. These results are always obtained 



from projected images. A plausible conversion to the 
3- dimensi onal fractal dimension, D ~ Dp + 1 is discussed 
bv lBeecbl (jT992i) . However, this relation holds rigorously 
only for the analysis of slices through an isotropic 
3-dimensional monofractal, while interstellar clouds 
could be multifra ctals (e.g.. IChappell fc Scald [2001). 
As discussed by IStutzki et al.l ( 19981 ) and shown by 
ISanchez et al.l ()2005[ ). Dp + 1 can be different from the 
3-dimensional fractal dimension for projected images. 

We applied the perimeter area method to projections 
along the x-, y- and z-axis, as well as to slices at a; = 0, 
y = and z = of our simulation data cubes. The results 
are presented in Figure [TT] for the projections (top panel) 
and the slices (bottom panel). Best-fit power laws to the 
projected data yield Dp ^1.36 for solenoidal forcing and 
Dp 1.18 for compressive forcing, whereas for the slices 
we find Dp^lA3 and Dpr^ 1.28, respectively. Thus, we 
find that the projections yield perimeter area dimensions 
systematically smaller than the slices for both forcings 
(Table 1). In agreement with the results obtained from 
the density Fourier spectra and the A-variance analysis, 
the perimeter area method yields a significantly smaller 
perimeter area dimension for compressive forcing com- 
pared to solenoidal forcing. 

4. CONCLUSIONS 

We have presented results of two high-resolution 
(1024^ grid cells) hydrodynamic simulations of super- 
sonic isothermal turbulence driven to rms Mach num- 
bers M ~ 5.5. The first simulation uses the typically 
adopted solenoidal (divergence-free) forcing to excite tur- 
bulent motions, whereas the second one uses compressive 
(curl- free) forcing. We have shown that compressive forc- 
ing yields much stronger density contrasts compared to 
solenoidal forcing for the same rms Mach number. This 
implies that the turbulence production mechanism leaves 
a strong imprint on compressible turbulence statistics, 
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especially altering the density statistics. Our results par- 
ticularly suggests that the mixture of solenoidal and com- 
pressive modes of the turbulence forcing must be taken 
into account. We summarize our results as follows: 

• The velocity Fourier spectra exhibit power laws 
in the inertial range for solenoidal and compres- 
sive forcing. The slopes obtained for both forcings 
are significantly steeper (~ 1.9) compared to the 
Kolmogorov slop e (5/3), in agreement with previ- 
ous s tudies fe.g.. ,Kritsuk et a l. 2007; Schmid t et al.l 
l2009l ) and in agreement with velocity dispersion 
to size relat ions infe rred from observations (e.g., 
LarsonlUQSUfF algaro ne et al.lll992l : iHever fc BruntI 
20041: iPadoan et al.ll20'06i r 

• From the integral of the density fluctuation Fourier 
spectra and from the asymptotic behavior of 
the 2nd order density structure function, we ob- 
tained the standard deviation of the density dis- 
tribution Up. Compressive forcing yields a stan- 
dard deviation ~ three times larger compared to 
solenoidal forcing, in agreement with the results 
found in our previous study a nalyzing density 
proba bility distribution functions (|Federrath et all 
20081) and in agreement with t h e studies by 
Passot fc Vazaucz-Semadenil (jl998f ). 'Kritsu k et al.l 
(|2007r) . iBectz et al, (,2008[ l and iSchmidt et al.l 

• The density fluctuation Fourier spectra are signifi- 
cantly steeper for compressive forcing in the inertial 
range compared to solenoidal forcing. Consistent 
results were obtained using complementary anal- 
ysis methods, i.e . , by comparing the A-variances 
(jOssenkopf et al.l l2008f ) and the 2nd order struc- 
ture functions of the density field. Our estimates of 
density spectra for solenoidal forcing are in agree- 
ment with previous studies, e.g., the weakly mag- 
netized super- Alfveni c su personic MHD mod els by 
IPadoan et all (|2004af ) and lKowal et al.l (|2007D . and 
consistent wi t h the hydrodynamic estimates by 
iKritsuk et all ()2006l ). Although a comparison with 
observational results must be regarded with cau- 
tion due to systematic uncertainties, our results for 
solenoidal and compressive forcing are in the range 
of inferred scaling exponents by observations (e.g., 
iBensch et al.ll200lD . 

• From the scaling of the density fluctuation Fourier 
spectra and the A-variance applied to the 3- 
dimensional data and applied to 2-dimcnsional pro- 
jections, we obtained f ractal Hurst expone nts fol- 
lowing the analysis bv IStutzki et al.1 (|1998D . This 
implies fractal box counting and fractal perimeter 
area dimensions significantly smaller for compres- 
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sive forcing compared to solenoidal forcing (see Ta- 
ble 1). 

• We analyzed the density structure us ing the fractal 
mass size method as introduced by IKritsuk et al.l 
(|2007f ). Compressive forcing yields a smaller frac- 
tal mass dimension compared to solenoidal forc- 
ing. The mass size method is, however, particu- 
larly sensitive to the temporal fluctuations of den- 
sity peaks. Given the large uncertainties, our 
results using this meth od are roughl y cons istent 
with the estimate s bv IKritsuk et all (|2007f) and 
iKowal fc LazarianI ()2007^ 

• We analyzed the fractal density structure using the 
box counting method described in Section 12.31 and 
the perimeter area method (Section 12. 5p applied 
to projections and slices. We recover the signifi- 
cant differences between solenoidal and compres- 
sive forcing inferred from the density spectra and 
A-variance analysis. However, the box counting 
dimension varies strongly with the defining den- 
sity threshold. The perimeter area dimensions ob- 
tained from slices are roughly consistent with the 
computed perimeter area dimensions from the A- 
variance given the systematic uncertainties (of or- 
der ~0.1 for fractal dimension estimates) compar- 
ing different methods. The range of fractal di- 
mensions obtained is consistent with the o b serva - 
tions analyzed by El megreen fc Falgarond (|1996l ) 
suggesting an overall fractal dimension of interstel- 
lar clouds in the range Z?'~2.3 ± 0.3. 
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TABLE 1 

Power-law exponents and fractal dimension estimates comparing solenoidal and compressive forcing 





symbol 


relation 


solenoidal forcing 


compressive forcing 


ID power-law index for 


6 




l.oD 


]^ g4±0.05 


ID power-law index for density 






0,7g±o.O6 


1 44±0.23 


3D power-law index for density 
fluctuation Fourier spectra 





P3D(fe) OC fc-" OC fe-<" + ^' 


2 78±o.06 


344±o.23 


derived Hurst exponent 

box counting dimension 

perimeter area dimension 


H 

Db 
Dp 


H = (/3-2)/2 
Db =3- H 
Dp = 2- H 


0.39±°"3 

2 gi±0.03 

1.61±°"^ 


Q 72±o i2 
2 28±o.i2 
1 28±oi2 


power-law index for A-variance 

applied to 3D data 







2 89±o.06 


344±o.i9 


derived Hurst exponent 

box counting dimension 

perimeter area dimension 


H 

Db 
Dp 


H = 0- 2)/2 
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Fig. 1. — Top panels: Column density computed along the 2;-axis in units of the mean column density for solenoidal forcing (left), and 
compressive forcing (right) at a randomly picked time t = 5T in the regime of statistically stationary compressible turbulence. Both maps 
show 4 orders of magnitude in column density with the same scaling for direct comparison of solenoidal and compressive forcing at rms 
Mach number ~5.5. Bottom panels: Same as top panels, but slices through the density field at z = 0. Compressive forcing yields stronger 
density enhancements and larger voids compared to solenoidal forcing. 
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t/T 

Fig. 2. — Bottom: rms Mach number JVi as function of the dynamical time T for 256^, 512'' and 1024^ numerical grid resolution. Top: 
Minimum and maximum density as function of the dynamical time T. At ^ 2T, a statistically stationary state was reached for both 
solenoidal (sol) and compressive (comp) forcing. Consequently, we can use all the available 81 snapshots within 2 < t/T < 10 for averaging 
statistical measures (e.g., Fourier spectra, structure functions. A- variance, fractal perimeter area, box counting and mass size analysis) to 
improve statistical significance and to compute corresponding Icr temporal fluctuations. Note that on average, the maximum density is 
almost ~ 10 times larger for compressive forcing compared to solenoidal forcing, alt hough the rms Mach number is roughly the same in both 
cases. The maximum density is subject to strongly intermittent fluctuations fe.g.. lKritsuk et al.ll2007h for both solenoidal and compressive 
forcing. 
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Fig. 3. — Kinetic energy Fourier spectra E{k) compensated by corresponding to Burgers turbulence (upper curves), and density 
fluctuation Fourier spectra compensated by k (lower curves) for solenoidal and compressive forcing respectively. Power-law fits in the 
inertial range 5 < A: < 15 are shown as thin solid lines. The velocity power spectra exhibit only marginal diff'erences between solenoidal and 
compressive forcing. The scaling of the density power spectra on the other hand differs significantly comparing both forcings. Accordingly, 
the stepper density power spectrum for compressive forcing leads to a smaller fractal box coverage dimension Dt, ~ 2.28 compared to 
Di,~2.61 for solenoidal forcing. 
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Fig. 4. — Top panel: Numerical resolution comparison of density fluctuation Fourier spectra for compressive forcing. At 512^ and 1024^, 
the spectra are almost converged with resolution, whereas the 256^ run deviates by ~ 15%. Bottom panel: Density fluctuation Fourier 
spectra at 1024'^ in comparison with its resampled version to 256^ cells. The resampling had virtually no influence on our results within 
the inertial range. 
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iogio(i/L) 

Fig. 5. — Top panel: A-variance analysis for the 3-dimensional data set resampled to a resolution of 2 56^ grid cells. Bottom panel: 
A-variance applied to the 2-dimensional projections of the 1024^ data set. As shown bv lStutzki et al] I I1998I ). the power-law scaling within 
the inertial range is preserved upon projection and agrees with the scaling of the Fourier power spectra of Figure|3]within the uncertainties 
from temporal fluctuations (see Tab. 1). 
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Fig. 6. — Second order structure functions of the density p for solenoidal and compressive forcing. The absolute values of the structure 
functions are in agreement with the measures of the power spectra and PDFs. The inertial range scaling, however, agrees with the Fourier 
spectra for compressive forcing only. Solenoidal forcing exhibits a density Fourier spectr um with powe r-law exponent a < 1 {Fig.|3]l, which 
precludes the determination of the power-law exponent via structure function analysis ijStutzki et al.. 1998. Appendix B). 
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Fig. 7. — Log-log plot of the mass M(l) compensated by obtained by the mass size method described in Section l2.4l for solenoidal and 
compressive forcing. As in Kowal & Lazarian (2007), a horizontal straight line therefore corresponds to a fractal mass dimension Dm = 3. 
Power-law fits in the inertial range yield a fractal mass dimension Dm ~ 2. 11 for solenoidal forcing and D,7i'^2.03 for compressive forcing. 
However, the difference between solenoidal and compressive forcing inferred from the mass size method is not significant due to the large 
uncertainties (error bars). These large uncertainties are a result of the strong temporal fluctuations of the maximum density (Fig. [2]l, since 
the mass size method relies on the density peaks, i.e., all cells with density p > pmax/2 are used as a basis for computing the mass size 
relation M{1). The dotted lines show M i"^ using Pmax only, i.e., without averaging over cells with p > Pmax/2. There is a systematic 
decrease of M(l) with decreasing averaging threshold, which does not affect the inertial range scaling as long as the averaging is performed 
over all cells with p > Pmax/2. 
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Fig. 8. — Log-log plot of -/V(Z) compensated with obtained by the box counting analysis described in Section 12.31 for solenoidal and 
compressive forcing respectively. A horizontal straight line would correspond to a box dimension of = 3. The power-law exponents D^^ 
obtained from fits within the inertial range are drawn as straight lines. Compressive forcing yields a significantly smaller box dimension 
Di,~2.51 compared to solenoidal forcing (Z)i,~2.68). 





Fig. 9. — Isosurface plots of the fractal structure defined by the mean density for solenoidal (ie/t) and compressive forcing (right). The pres- 
ence of hierarchical visual complexity indicates a fractal structure that can be analyzed with the box counting method (.Mandelbrot fc Frarn^ 
[200^ . 
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Fig. 10. — Shows the dependence of the box counting dimension on the threshold density ptji defining the fractal structure. Computing 
the fractal box counting dimension for t = logj^Q(pth/po) = ~li Oi li 2 reveals the strong dependence of the inferred fractal dimension on 
the defining density threshold. Left panel: solenoidal forcing; Right panel: compressive forcing. 
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Fig. 11. — Top panel: Perimeter area method applied to projections along x, y and z for solenoidal and compressive forcing respectively. 
Bottom panel: Same as top panel but for slices at a; = 0, y = and 2 = 0. The perimeter is given in units of the numerical cell size 
A = L/1024. 



